library(car)  
library(plyr) 
library(dplyr) 
library(foreign) 
library(ggplot2) 
library(lme4) 
library(mapproj) 
library(maps)
library(pbapply) 
library(reshape2)
library(rms)
library(stringr) 
library(survey) 
library(devtools) 
library(tidyr)

FileName <- function (path=NULL, name=NULL, ext=NULL, replace=FALSE) {
  p <- paste(path, collapse = "")
  n <- paste(name, collapse = "")
  e <- paste(".", ext, sep = "")
  d <- format(Sys.Date(), "%y%m%d") 
  for (l in seq_along(letters)) {
    if (l > 1) old_file_name <- file_name
    file_name <- paste0(p, d, paste(n, letters[l], sep="-"), e)
    if (!any(grepl(file_name, list.files()))) {
      if (replace) file_name <- ifelse(l > 1, old_file_name, file_name)
      break
    }
  }
  cat("\nFile name:", file_name, "\n\n")
  return(file_name)
}

mypdf <- function(name, ...) {
  date <- format(Sys.Date(), "%y%m%d")
  pdf(paste(paste(name, collapse=""), date, ".pdf", sep=""), ...)
}
mysw <- function (...) {
  setwd(paste(..., sep = "/"))
}
# Multiple plot function
multiplot <- function(..., plotlist = NULL, file, cols = 1, layout = NULL) {
  require(grid)

                                        # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

                                        # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
                                        # Make the panel
                                        # ncol: Number of columns of plots
                                        # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), 
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }

  if (numPlots == 1) {
    print(plots[[1]])

  } else {
                                        # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

                                        # Make each plot, in the correct location
    for (i in 1:numPlots) {
                                        # Get the i, j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout ==  i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, 
                            layout.pos.col = matchidx$col))
    }
  }
}


############
## Figure 2: Map of climate concern
############


### PLOTS ###
theme_clean <- function (base_size = 12) { 
  require(grid) # Needed for unit() function
  theme_grey(base_size) %+replace%
  theme(
      axis.title = element_blank(), 
      axis.text = element_blank(), panel.background = element_blank(), 
      panel.grid = element_blank(), axis.ticks.length = unit(0, "cm"), 
      axis.ticks.margin = unit(0, "cm"), 
      panel.margin = unit(0, "lines"), 
      plot.margin = unit(c(0, 0, 0, 0), "lines"), complete = TRUE
  ) 
}


POtoCode <- function (pos, dta = st.info) {
  return(as.integer(dta$ICPSRCode[match(pos, dta$POAbrv)]))
}
CodeToName <- function (stcodes, dta = st.info) {
  return(as.character(dta$Name[match(stcodes, dta$ICPSRCode)]))
}
FileName <- function (path = NULL, name = NULL, ext = NULL) {
  p <- paste(path, collapse = "")
  n <- paste(name, collapse = "")
  e <- paste(".", ext, sep = "")
  d <- format(Sys.Date(), "%y%m%d") 
  file_name <- paste(p, d, n, "-a", e, sep = "")  
  if (any(grepl(file_name, list.files()))) {
    for (ll in letters[-1]) {
      file_name <- paste(p, d, paste(n, ll, sep = "-"), e, sep = "") 
      if (!any(grepl(file_name, list.files()))) break
    }
  }
  cat("\nFile name:", file_name, "\n\n")
  return(file_name)
}


states <- read.csv(file = "states_icpsr.csv")

states$NAME <- gsub(" ", "", states$NAME)
states$NAME[states$NAME == "DistrictofColumbia"] <- "District of Columbia"
states$NAME[states$NAME == "NewHampshire"] <- "New Hampshire"
states$NAME[states$NAME == "NewJersey"] <- "New Jersey"
states$NAME[states$NAME == "NewMexico"] <- "New Mexico"
states$NAME[states$NAME == "NewYork"] <- "New York"
states$NAME[states$NAME == "NorthCarolina"] <- "North Carolina"
states$NAME[states$NAME == "NorthDakota"] <- "North Dakota"
states$NAME[states$NAME == "RhodeIsland"] <- "Rhode Island"
states$NAME[states$NAME == "SouthCarolina"] <- "South Carolina"
states$NAME[states$NAME == "SouthDakota"] <- "South Dakota"

states$NAME[states$NAME == "WestVirginia"] <- "West Virginia"

states <- states[, c(1, 2, 6)]
states$state <- tolower(states$NAME)

climate_scale2<- read.dta(file="climate_analysis_states.dta")
climate_scale3<-climate_scale2

climate_scale3<-dplyr::group_by(climate_scale3, year, abb) %>%
	dplyr::summarise_each(funs(mean) )
	
opinion_quantiles<-dplyr::group_by(climate_scale2, year, abb) %>%
	summarise(
		mass_climate_concern_lower=quantile(mass_climate_concern, probs=0.05),
		mass_climate_concern_median=quantile(mass_climate_concern, probs=0.5),
		mass_climate_concern_upper=quantile(mass_climate_concern, probs=0.95)
	)
climate_scale3<-merge(climate_scale3, opinion_quantiles, by=c("year", "abb"))	

##############
### Figure 3: State map of climate concern
##############

climate_scale <- climate_scale3
climate.dta <- climate_scale
year.median <- tapply(climate.dta$mass_climate_concern_median, climate.dta$year, median, na.rm = T)
liberalism.map <- map_data("state")

climate.dta <- climate.dta[climate.dta$year ==
                                             2001  |
                                         climate.dta$year ==
                                             2009 |
                                         climate.dta$year ==
                                             2017,]


climate.dta$Year.Median[climate.dta$year == 2001] <-
  mean(climate_scale$mass_climate_concern_median[climate_scale$year == 2001])

climate.dta$Year.Median[climate.dta$year == 2009] <-
  mean(climate_scale$mass_climate_concern_median[climate_scale$year == 2009])

climate.dta$Year.Median[climate.dta$year == 2017] <-
  mean(climate_scale$mass_climate_concern_median[climate_scale$year == 2017])


climate.dta$demeaned <- climate.dta$mass_climate_concern_median-climate.dta$Year.Median

climate.dta <- mutate(climate.dta, StPO = abb, State = abb, 
                         state = tolower(abb))
climate.dta <- climate.dta[!is.na(climate.dta$year), ]
climate.dta <- ddply(as.data.frame(climate.dta), ~year, mutate,
                        liberalism0 = scale(median)[, 1])

years <- c( 2001, 2009, 2017)

for (i in 1:3){
  climate.dta$policy0[climate.dta$year == years[i]] <- 
    scale(climate.dta$mass_climate_concern_median[climate.dta$year == 	
                                    years[i]], center = T, scale = T)
}
liberalism.map <- merge(liberalism.map, states, by.x = "region", by.y = "state")
liberalism.map <- merge(liberalism.map, climate.dta, by.x = "STATEAB", by.y = "abb")
liberalism.map <- arrange(liberalism.map, group, order)



pdf("Figure3_map_climate_opinion.pdf", height = 2)
(p1 <- ggplot(data = liberalism.map, 
              aes(x = long, y = lat, group = group, fill = policy0))
 + geom_polygon(colour = "black", size = rel(.1))
 + coord_map("polyconic")
 +  viridis::scale_fill_viridis(direction=-1)
 #+ scale_fill_gradient2(name = "Liberalism", 
 #                       low = 'red', mid = "white",high = 'blue')

 + facet_wrap(~year, ncol = 3)
 + theme_clean()+guides(fill = FALSE)
 )
dev.off()


##############
### Figure C1: Validation Analysis using Yale Estimates
##############

library(dplyr)

YaleEstimates<-read.csv(file="nclimate2583_state.csv")
YaleEstimates$Estimate_year<-2012
    
climate_scale3<-merge(climate_scale3, YaleEstimates,by.x=c("abb","year"), by.y=c("State_abb","Estimate_year"), all.x=T)
 cor(climate_scale3$mass_climate_concern_median, climate_scale3$x65_happening, use="complete.obs")
  cor(climate_scale3$mass_climate_concern_median, climate_scale3$x67_human, use="complete.obs")
 cor(climate_scale3$mass_climate_concern_median, climate_scale3$x73_consensus, use="complete.obs")
 cor(climate_scale3$mass_climate_concern_median, climate_scale3$x78_worried, use="complete.obs")
 cor(climate_scale3$mass_climate_concern_median, climate_scale3$x82_harmUS, use="complete.obs")
climate_scale3$constant<-1

climate_scale_long <- gather(climate_scale3, yale_question, yale_mrp, x65_happening:x82_harmUS, factor_key=TRUE)

climate_scale_long$yale_label<-NA
climate_scale_long$yale_label[climate_scale_long$yale_question=="x65_happening"]<-"Global Warming is Happening"
climate_scale_long$yale_label[climate_scale_long$yale_question=="x67_human"]<-"Climate Change is Human-Caused"
climate_scale_long$yale_label[climate_scale_long$yale_question=="x73_consensus"]<-"Believe Most Scientists think Global Warming Happening"
climate_scale_long$yale_label[climate_scale_long$yale_question=="x78_worried"]<-"Worried about Climate Change"
climate_scale_long$yale_label[climate_scale_long$yale_question=="x82_harmUS"]<-"Global Warming Will Harm United States"
 
climate_scale_long_grouped = group_by(climate_scale_long, yale_label)
cors<-dplyr::summarize(climate_scale_long_grouped, cor(mass_climate_concern_median, yale_mrp,use = "complete.obs"))
colnames(cors)[2]<-"cor"

validation1 <- (ggplot(data = subset(climate_scale_long, year==2012), aes(x = mass_climate_concern_median, y =
                                         yale_mrp, label=abb))
      + theme_bw()
      + theme(axis.title.y = element_text(size = 16))
     + geom_text(size=3)
        + geom_smooth(method = "lm", se = FALSE)
    + facet_wrap(~yale_label, ncol = 2)
      + ylab("Howe et al (2015) MRP Estimates of State-Level Opinion about Climate Change")
      + xlab("Climate Concern Index (2012)")
      + geom_text(data = cors, aes(label = paste("r = ", round(cor, 2), sep =
                                                     "")), size = 6, x
                  = 1.25, y =35)
      +theme(plot.title = element_text(size = 20), axis.title.y =
                     element_text(size = 8))
     )


pdf("FigureC1_validation_YaleStudy.pdf", height = 8)
    validation1
    dev.off()



